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Abstract. We use large-scale three-dimensional simulations of supersonic Euler turbulence to study 
the physics of a highly compressible cascade. Our numerical experiments describe non-magnetized 
driven turbulent flows with an isothermal equation of state and an rms Mach number of 6. We 
find that the inertial range velocity scaling deviates strongly from the incompressible Kolmogorov 
laws. We propose an extension of Kolmogorov's K41 phenomenology that takes into account 
compressibility by mixing the velocity and density statistics and preserves the K41 scaling of the 
density-weighted velocity v = p '/^m. We show that low-order statistics of v are invariant with respect 
to changes in the Mach number. For instance, at Mach 6 the slope of the power spectrum of v 
is —1.69 and the third-order structure function of v scales linearly with separation. We directly 
measure the mass dimension of the "fractal" density distribution in the inertial subrange, D„, « 2.4, 
which is similar to the observed fractal dimension of molecular clouds and agrees well with the 
cascade phenomenology. 
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INTRODUCTION 

In the late 1930's, Kolmogorov clearly realized that chances to develop a closed purely 
mathematical theory of turbulence are extremely low [1].^ Therefore, the basic approach 
in [2, 3] (usually referred to as the K41 theory) was to rely on physical intuition and 
formulate two general statistical hypotheses which describe the universal equilibrium 
regime of small-scale fluctuations in arbitrary turbulent flow at high Reynolds number. 
Following the Landau (1944) remark on the lack of universality in turbulent flows [4], 
and with information extracted from new experimental data, the original similarity hy- 
potheses of K41 were then revisited and refined to account for intermittency effects 
[5, 6, 7]. While the K41 phenomenology became the cornerstone for all subsequent de- 
velopments in incompressible turbulence research [e.g., 8], there was no similar result 
established for compressible flows yet [9, 10]. Historically, compressible turbulence re- 
search, preoccupied with a variety of specific engineering applications, was generally 



"An understanding of solutions to the [incompressible] Navier-Stokes equations" yet remains one of the 
six unsolved grand challenge problems nominated by the Clay Mathematics Institute in 2000 for a $1M 
Millennium Prize [http://www.claymath.org/millennium/]. 




0.5 1 1.5 2 ' 1.4 1.6 1.8 2 2.2 2.4 



FIGURE 1. Time average compensated power spectra (left) and third-order transverse structure func- 
tions (right) for velocity u and mass- weighted velocities v = p^^^u and w = p^^^u. The statistics of v 
clearly demonstrate a K41-like scaUng. Notice strong bottieneck contamination in the spectra at high 
wavenumbers. 



lagging behind the incompressible developments.^ The two major reasons for this time 
lag were an additional complexity of analytical treatment of compressible flows and a 
shortage in experimental data for super- and hypersonic turbulence. In this respect, al- 
though limited to relatively low Reynolds numbers, direct numerical simulations (DNS) 
of turbulence (pioneered by Orszag and Patterson [14]) have occupied the niche of ex- 
periments at least for the most simple flows. One particularly important advantage of 
DNS is an easy access to variables that are otherwise difficult to measure in the labora- 
tory or treat analytically. 

A traditionally straightforward approach to data analysis from DNS of compressible 
turbulence includes computation of the "standard" statistics of velocity fluctuations. In 
addition, the diagnostics for density fluctuations are also computed and discussed as the 
direct measures of compressibility. Quite naturally, both density and velocity statistics 
demonstrate strong dependence on the Mach number ^ in supersonic (^ E [1,3]) and 
hypersonic (^ > 3) regimes, while the variations in turbulent diagnostics at sub- or 
transonic Mach numbers are rather small. For instance, at ^ 1 the velocity power 
spectrum closely follows the K41 scaling and the third order velocity structure functions 
scale roughly linearly with separation [15]. The density power spectrum in weakly 
compressible isothermal flows scales as ~ k~'^/^ [16], at ^ a; 1 it scales as ~ k~^'^ 
[20], and at ^ 6 the slope is -1.07 [20]. 

Based on the data from numerical experiments, it is well established that: (i) the 
velocity power spectra tend to get steeper as the Mach number increases, reaching the 
Burgers slope of —2 asymptotically [18, and references therein]; (ii) the density power 
spectra instead get shallower at high Mach numbers, approaching a slope of —1 or even 
shallower [17]; (iii) the density PDF in isothermal turbulent flows is well represented by 
a lognormal distribution [18, and references therein]; (iv) the dimensionality of the most 



A reasonable measure of the delay is 60+ years passed between the appearance of incompressible 
Reynolds averaging [11] and mass-weighted Favre averaging for fluid flows with variable density [12], 
although see [13] for references to a few earlier papers that dealt with density-weighted averaging. 



singular velocity structures increases from Ds ^ ~ 1 in a subsonic regime to ~ 2 
in highly supersonic [19]; (v) the mass dimension of the turbulent structures decreases 
from Dm = 3 in weakly compressible flows to ~ 2.5 in highly compressible [20]. 

How can we combine these seemingly disconnected pieces of information into a 
coherent physical picture to improve our understanding of compressible turbulence? 
One way to do this is to consider a phenomenological concept of a lossy compressible 
turbulent cascade that would asymptotically match the incompressible Kolmogorov- 
Richardson energy cascade [2, 21] in the limit of very low Mach numbers. Since 
incompressible turbulence represents a degenerate case where the density is uncorrelated 
with the velocity, the phenomenology of the compressible cascade must include this 
correlation. This essentially means that instead of velocity u, which is a single key 
ingredient of the K41 laws, one needs to consider a set of mixed variables, p^/^u, where 
p is the density and rj can take values 1, 2, or 3, depending on the statistical measure of 
interest [20]. For instance, if one is studying the scale-by-scale kinetic energy budget in 
a compressible turbulent flow, a mixed variable power spectrum with t] = 2 would be an 
appropriate choice. To deal with the kinetic energy flux through the hierarchy of scales 
within the inertial range, the key mixed variable would be the one with 7] = 3. 

How will these mixed statistics scale in the inertial range of highly compressible tur- 
bulent flows? Will their scaling depend on the Mach number? Can the K41 phenomenol- 
ogy be extended to cover hypersonic turbulent flows? These and other related questions 
are in detail discussed in [20] based on Euler simulations of driven isotropic supersonic 
turbulence with the Piecewise Parabolic Method [22] and with resolution up to 2048-^ 
grid points. In this paper we present the highlights of the compressible cascade phe- 
nomenology verified in [20]. 

SCALING, STRUCTURES, AND INTERMITTENCY 

Nonlinear interactions transfer kinetic energy supplied to the system at large scales 
through the inertial range with little dissipation. Let us assume that the mean volume 
energy transfer rate in a compressible fluid, pu^u/l, is constant in a statistical steady 
state [e.g., 23]. If this is true, then 

vP = {p^i^u)p ^^p|^ (1) 

for an arbitrary power p and, with the standard assumption of self- similarity of the 
cascade, the structure functions (SFs) of mixed variable v for compressible flows should 
scale in the inertial range as 

= (\v{r + E)-v{r)\P)r^EPl\ (2) 

In the limit of weak compressibility, the scaling laws (2) will reduce to the K41 
results for the velocity structure functions. The scaling laws S^p{t) ^ I^p, where C,p = 
p/3 are not necessarily exact. As the incompressible K41 scaling, they are subject to 
"intermittency corrections", e.g. ll,p = p/'i-\- [5]- The only exception is, perhaps, the 
third order relation for the longitudinal velocity SFs, which is exact in the incompressible 
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FIGURE 2. Gas mass M(^) as a function of the box size t (left). The mass dimension is defined as 
the log-log slope of M{£), see eq. (4). Relative exponents for structure functions of the transverse modified 
velocities v versus order p and two hierarchical structure models with different parameters [HS 1 & HS2, 
6] that fit the data for p G [0, 3] (right). Also shown are model predictions for the Kolmogorov-Richardson 
cascade [K41, 2, 3], for intermittent incompressible turbulence [SL94, 6], for "burgulence" [Burg, 24], 
and for the velocity fluctuations in supersonic turbulence [B02, 25]. 



case and is known as the four-fifth law [3]. Our focus here is mostly on the low order 
statistics (p < 3) for which the corrections are small. Since the power spectrum slope is 
related to the exponent of the second order structure function, the K41 slope of 5/3 is 
expected to hold for v = p ^^^u in the compressible case. 

Figure 1 shows the power spectra of u, v, and w = p^^^u and the corresponding third- 
order transverse structure functions based on the simulations at Mach 6 [26, 20]. The 
power spectrum E(^) and the structure function of v clearly follow the K41 scaling: 
L ~ k~^-^^ and ^3 ~ £^-^^ [20], while the velocity power spectrum (^{k) and struc- 
ture function have substantially steeper-then-K41 slopes: —1.95 and 1.29 [26]. At the 
same time, the kinetic energy spectrum E ~ k-^-^^ is shallow and both solenoidal and 
dilatational components of w have the same slope implying a single compressible en- 
ergy cascade with strong interaction between the two components [20]. These results 
based on the high dynamic range simulations lend strong support to the scaling relations 
described by eq. (2) and to the conjecture from which they were inferred. Previous sim- 
ulations at lower resolution did not allow to measure the absolute exponents reliably due 
to insufficient dynamic range and due to the bottleneck contamination [27]. 

In 195 1, von Weiszaker [28] introduced a phenomenological model for scale-invariant 
hierarchy of density fluctuations in compressible turbulence described by a simple 
equation that relates the mass density at two successive levels to the corresponding scales 
through a universal measure of the degree of compression, a, 

p„/p„_i = (4/4-i)~^". (3) 

The geometric factor a takes the value of 1 in a special case of isotropic compression 
in three dimensions, 1/3 for a perfect one-dimensional compression, and zero in the 
incompressible limit. From equations (1) and (3), assuming mass conservation. Fleck 
[29] derived a set of scaling relations for the velocity, specific kinetic energy, density. 




FIGURE 3. Coherent structures in Mach 6 turbulence at resolution of 1024^. Projections along the 
minor axis of a subvolume of 700 x 500 x 250 zones for the density (upper left), the enstrophy (upper 
right), the dissipation rate (lower left), and the dilatation (lower right). The logarithmic grey-scale ramp 
shows the lower values as dark in all cases except for the density. The inertial subrange structures 
correspond to scales between 40 and 250 zones and represent a fractal with ~ 2.4. The dominant 
structures in the dissipation range (i < 30A) are shocks with = 2. [Reprinted from [20].] 



and mass: 

where all the exponents depend on the compression measure a which is in turn a 
function of the rms Mach number of the turbulent flow. We can now use the data from 
numerical experiments to verify the scaling relations (4). Since the first-order velocity 
structure function scales as £^-^^ [20], we can estimate a for the IVIach 6 flow, a ^ 0.21. 
Using the last relation in (4), we can calculate the mass dimension for the density 
distribution, ~ 2.38. It is indeed consistent with our direct measurement of the mass 
dimension for the same range of scales. Dm ~ 2.39, see Fig. 2. 

In strongly compressible turbulence at IVIach 6, the density contrast between superson- 
ically moving blobs and their more diffuse environment can be as high as 10^. The most 
common structural elements in such highly fragmented flows are nested bow-shocks 
[17]. Figure 3 shows an extreme example of structures formed by a collision of counter- 
propagating supersonic flows. On small scales within the dissipation range, these struc- 
tures are characterized by = 2, while within the inertial range ~ 2.4 (Fig. 2, left). 



The hierarchical structure (HS) model 



Cp/C3 = rp+c{i-pn (5) 

[6] provides good fits to the data for the mass-weighted velocity v (see Fig. 2, right). 
Here the codimension of the support of the most singular dissipative structures 

C=3-D,,,= (l-37)/(l-j83). (6) 

If the fit is limited to p e [0, 3], two sets of model parameters j8 and 7 are formally 
acceptable (models HSl and HS2 in Fig. 2). The best-fit parameters of the HSl model: 
j8f = 1/3 (a measure of intermittency), = (a measure of singularity of structures), 
and Ci = 1 .5 correspond to a hybrid between the B02 model for the velocity fluctuations 
(/3J02 = 1/3, Yboi = 1/9) [25] and the Burgers' model (^surg = 0, yeurg = 0) [24]. The 
HS2 model (jSj = 1/6, 72 = 1/9, and C2 = 0.8) provides a fit of roughly the same 
quality for p G [0, 3], but overestimates the scaling exponents l^p at p > 4. Since the 
level of uncertainty in the high order statistics remains high even at a resolution of 1024^ 
grid points, larger dynamic range simulations are needed to distinguish between the two 
options. 

If the HSl option is confirmed, then Mach 6 turbulence is more intermittent than 
incompressible turbulence (jSj' < /3|^94 = 2/3) and has the same degree of singularity 
of structures as burgulence. The singular dissipative structures with fractal dimension 
Dg^y = 1.5 can be conceived as perforated sheets reminiscent of the Sierpinski sieve. If 
the HS2 option is justified, then turbulence at Mach 6 is even more intermittent, but 
has the same degree of singularity of structures as incompressible turbulence (Y2. = 
YsL94-)- In this case the fractal dimension of the most singular structures, Ds^v = 2.2, 
is slightly higher than in the B02 model [25]. Formally, it is also possible that both 
types of structures are present in highly compressible turbulence, implying multiple 
modulation defects and a compound nature of Poisson statistic [cf. 30]. In this case, 
a linear combination of HS 1 and HS2 models would describe the high order exponents 
best. 

CONCLUSION 

Using large-scale Euler simulations of supersonic turbulence at Mach 6 we have demon- 
strated that there exists an analogue of the K41 scaling laws valid for both weakly and 
highly compressible flows. The mass-weighted velocity v = p^^^u - the primary vari- 
able governing the energy transfer through the cascade - should replace the velocity u 
in intermittency models for compressible flows at high Mach numbers. 
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